Generalized London free energy for high- Tc vortex lattices 
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We generalize the London free energy to include four-fold anisotropies which could arise from d- 
wave pairing or other sources in a tetragonal material. We use this simple model to study vortex 
lattice structure and discuss neutron scattering, STM, Bitter decoration and jiSK experiments. 



The London free energy provides a very simple way 
of studying the vortex lattice in an extreme type II su- 
perconductor. The conventional isotropic model pre- 
dicts a hexagonal vortex lattice with an arbitrary ori- 
entation relative to the ionic lattice. Recent neutron 
scattering Q and STM experiments on high-Tc com- 
pound YBa2Cu307_5 (YBCO) revealed vortex lattices 
with centered rectangular symmetry and a specific ori- 
entation with respect to the ionic lattice. It has been 
proposed that this effect can arise from additional quar- 
tic derivative terms in the Ginzburg-Landau (G-L) free 
energy or, alternatively, from including two or more 
order parameters (such as d and s) in the G-L free energy 
with derivative mixing terms reflecting the ionic lattice 
symmetry . Such models predict interesting effects 

in the behavior of the various order parameters in the 
vortex lattice. However, these models contain a large 
number of unknown parameters and are rather cumber- 
some to work with numerically. Another approach [^jl3| 
to the macroscopic effects of c?-wave pairing takes into ac- 
count the generation of quasiparticles near the gap nodes 
due to current flow and thermal excitation. This leads to 
a non-linear relationship between supercurrent and su- 
pcrfluid velocity which becomes singular at T ^ 0. 

The purpose of this letter is to present a simple and 
general approach to these effects based on a generaliza- 
tion of the London free energy to include anisotropy of 
four-fold symmetry, characteristic of a tetragonal ionic 
lattice. The number of new parameters is far smaller 
than in the G-L approach (a reasonable model contains 
only one new parameter) and numerical simulations are 
considerably easier. It provides a useful model to study 
vortex lattice structure, pinning by twin boundaries and 
the magnetic field distribution measured in /iSR experi- 
ments. The model is suitable to study the intermediate 
field region Hd <C -ff <C Hc2 which is experimentally 
most relevant but traditionally difficult to handle within 
the G-L theory. Furthermore, this approach can be ex- 
tended to r = where G-L theory breaks down and the 
supercurrent becomes singular. 

We now present a derivation of the generalized London 
model, starting from a G-L free energy density with both 
d and s order parameters |^,|l^,|l^ : 

/ - + arflrfp + 7,|ns|2 + ^a\lid\^ + + hySir 



+ 7,[(nys)*(nyrf) - (n,s)*(n,d) -f c.c.]. (i) 

Here H = — iV — e*A/hc and fi contains the quartic 
terms. We shall consider a case of a d-wave superconduc- 
tor in which s identically vanishes in zero magnetic field. 
In finite field (H > Hd) a small s-component with a 
highly anisotropic spatial distribution is nucleated in the 
vicinity of a vortex givingrise to non-triangular equilib- 
rium lattice structures [^|o) . Our strategy will be to sim- 
plify free energy (|l|) by integrating out this s-component 
in favor of higher order derivative terms in d. In this 
process some short length-scale information on the order 
parameter is lost but the magnetic field distribution is 
described accurately. Using its Euler-Lagrange equation 
s can be expressed to the leading order in (1 — T/Tc) as 

s={j,/as){Ill-Ul)d. (2) 

Substituting this into / gives the leading derivative terms 
in d of the form: 

/ = 7rf[|nrfP - (7.V7rf«s)|(n2 - H2)dn + ... (3) 

Various additional corrections to the free energy are ob- 
tained from integrating out s more accurately, taking into 
account the 7s|Hsp term and quartic terms. However 
these all involve higher powers of H or other terms that 
will not concern us. The coefficient of the second term 
has dimensions of (length)^; we will write it in the form 
where e = S^adjy / ctsJd) is a dimensionless param- 
eter which controls the strength of the s-d coupling and 
^ = \/jd/\ad\ is the G-L coherence length. We hence- 
forth assume e <g; 1. As we remark below, neutron scat- 
tering and STM experiments probably support this as- 
sumption. We note that a term of the form [(H^ — H^)(ip 
could arise without invoking s-d mixing from a systematic 
derivation of higher order terms in the G-L free energy 
starting with a BCS-like model and taking into account 
the square symmetry of the Fermi surface fl-|^,[l6|] . 

The free energy of Eq. ^ is not bounded below, ex- 
hibiting runaway behavior for rapidly varying d-fields. 
This is in fact cured by keeping additional higher deriva- 
tive terms that also arise from integrating out s. Sta- 
bility occurs for 7^ < 7^7^- In fact, the approximation 
of Eq. (^) will be sufficient for our purposes, yielding a 
local minimum which we expect would become a global 
minimum upon including the additional terms. 
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We now assume that the penetration depth A ^ ^. We 
may then assume that |d(r)| « do, the zero field equi- 
hbrium value, almost everywhere in the vortex lattice, 
except within a distance of O(^) of the cores. This gives 
the London free energy, 

^{dyVy-d.,v^f\}, (4) 

written in terms of the superfluid velocity, 

= V6' - (eV^c)!, (5) 

where 9 is the phase of d. 

The corresponding London equation, obtained by vary- 
ing with respect to A, is: 



- [ydy - xdx){dyVy - d,j,v^)W. 



(6) 



For many purposes it is very convenient to express v in 
terms of B and its derivatives, and then substitute this 
expression for v back into /i, giving an explicit expres- 
sion for /l as a functional of B only. For e = this 
gives 

t)<o) V X B/Bo, (7) 

where Bq = (/)o/27rA^ is of order Hd {(po = 2TThc/e* is 
the flux quantum) and 

/° = (l/87r)[S2 + A2(VxB)2]. (8) 

Here the penetration depth, for e = is Aq^ = 
8n^d{e*do/hc)^ . It is presumably not possible to solve 
Eq. in closed form for u as a function of i? for e 7^ 0. 
However, this can be done readily in a perturbative ex- 
pansion in e. The first order correction is: 



(9) 



with t/'^"^ given by Eq. (^. 
density, up to 0(e) is then: 

fL = fL + '-^[^{d.dyBr- 



The London free energy 



[dyBff/Bl]. 

(10) 



Note that we could have arrived at a similar conclusion 
by simply writing down all terms allowed by symmetry 
in /i, expanding in number of derivatives and powers 
of B. Square anisotropy is first possible in the fourth 
derivative terms. In principle, we should also include all 
isotropic terms to order B^ and V*. However, assuming 
that these have small coefficients, they will not be impor- 
tant. This result can also be obtained from considering 
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FIG. 1. a) Distribution of magnetic field in a vortex lattice 
for e = 0.3 and H = 6.8T, leading to an angle of /? ~ 74° . We 
use Ao = 1400A and k = Ao/^ = 68. b) Corresponding /iSR 
lineshape. 

generation of quasi-particles near gap nodes in a d-wave 
superconductor ||l^ , in a range of temperature and field 
where the supercurrent can be Taylor expanded in the 
superfiuid velocity |l7|. More generally, the quadratic 
and quartic terms in ([10[) have independent coefficients. 

The corresponding London equation is obtained by 
varying Jl with respect to B{r). For B along the z di- 
rection one obtains 



[1 - A^V^ + Ae\ie{d.cdyY]B - eQ[B] - 0, 



(11) 



where 



Q[B] = 2\ieB^\dl - dl)B + d^Bd, 
X [{d^Bf - [dyBf] 



dyBdy] 



(12) 



is the non-linear term arising from the last term in Eq. 
®- 

To get a feeling for the effect of the extra terms, con- 
sider a weak field which depends only on x or else only 
on {x + y). The solution of the linearized London equa- 
tion [(pd]) without the last term] gives an exponentially 
decaying field with A = Aq for variation along the x-axis 
but: 



A = Ao 




1/2 



(13) 
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for variation at 45° to the crystal axis ||l8[ . The penetra- 
tion depth is longer along the crystal axis. 

To determine vortex lattice structure we insert source 
terms p{r — rj) at the vortex core positions, Vj, on the 



right hand side of Eq. (11). The source terms reflect the 
topological winding of the phase angle and the reduction 
of the order parameter in the core W . A commonly used 
phenomenological form is [p9[: 



p{r) = (0o/2<2)g 



(14) 



It is straightforward to solve these equations numerically 
for the vortex lattice by an iterative method. We find 
that the quartic term makes a negligible contribution. 
(Contrary to naive expectation, it doesn't become more 
important with increasing applied field because the field 
becomes nearly constant in the vortex lattice when the 
applied field is large.) Thus to an excellent approxima- 
tion one may neglect Q[B\ in the London equation (|lT 
and the magnetic field may be written explicitly as: 



(15) 



Here the sum is over all wave-vectors in the reciprocal 
lattice and B is the average field. The lattice constant is 
determined by the condition that B^l = (pQ where n is the 
area of the unit cell. The lattice symmetry is determined 
by minimizing the Gibbs free energy gL — /l — HE /At:. 
We find that the flux lattice has centered rectangular 
symmetry, with principal axes aligned with the ionic crys- 
tal lattice, with an angle (3 between unit vectors which 
depends on e and the magnetic field. An example of 
such a centered rectangular lattice is shown in Fig. ^(a). 
In agreement with earlier results within G-L llfl] and 
Eilenberger |pT[| formalisms individual vortices are elon- 
gated along the crystalline axes. Figure |2|(a) shows the 
dependence of Gibbs free energy on (3 for various val- 
ues of e at fixed apphed field H = 400Bo ~ 6.8T. For 
e = minimum occurs for /3min = 60°, corresponding 
to a hexagonal lattice. As e increases /3min continuously 
increases and for sufficiently large e, the fiux lattice be- 
comes tetragonal with /3min = 90°. For (iyim 7^ 90°, 
there are always two solutions, related by a 90° rota- 
tion, in which the long axis of the centered rectangle is 
aligned with either the x ot y axis. The degeneracy is 
much larger for e = 0, when the flux lattice may have an 
arbitrary orientation relative to the ionic crystal lattice. 

Dependence of /3min on the applied field for various val- 
ues of e is displayed in Fig. |2|(b) . Clearly the anisotropic 
term becomes more important at larger fields. Our per- 
turbative elimination of v in favor of B breaks down when 
e and H are sufficiently large that /3min differs signifi- 
cantly from 60°. Furthermore, we might expect higher 
order corrections to (^) to be important in this regime. 
By fitting Fig. ^(b) to experimental data on tetragonal 
materials such as Tl2Ba2Cu06+(i (once such data become 
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FIG. 2. a) Gibbs free energy as a function of /3 for the 
same parameters as Fig. |l| and various values of e. Arrows 
indicate positions /3min of the minima and Go = -HVStt. b) 
Equilibrium angle /3min as a function of H for several values 
of e. 



available) one can directly assess the magnitude of e, the 
only unknown parameter in the model. 

Our analysis can be easily extended to take into ac- 
count effective mass (i.e. penetration depth) anisotropy. 
In a simple one-component G-L model, the derivative 
term is generalized to: 



(16) 



We restrict our attention to fields along the z-axis. Then 
the anisotropy can be removed by a rescaling of the x- 
coordinate and a corresponding rescaling of the mag- 
netic field. The coherence length and penetration depth 
anisotropics are the same: i,y/S,x = ^x/^y We will make 
the simplifying assumption that the higher derivative and 
mixed derivative terms in / are also simply modified by 
a rescaling by a common factor. It then follows that the 
ffux lattice shape is obtained by stretching along the x- 
axis by the factor X^/Xy We now obtain two possible 
vortex lattices, both of centered rectangular symmetry, 
aligned with the ionic lattice, with different angles, (3. 
(Relaxing our simplifying assumption may split the de- 
generacy between these two lattices.) On the other hand, 
when e = 0, we may rotate the hexagonal lattice by an 
arbitrary angle before stretching. This gives an infinite 
set of oblique lattices with arbitrary orientation. 

To compare theory with YBCO we should take into 
account twin boundaries which may also tend to align 
the vortex lattice by pinning vortices to the twin bound- 
aries, at ±45*^ to the x-axis. This effect competes with 
alignment to the ionic lattice which we have been dis- 
cussing. Only in the special case of a square vortex lat- 
tice does a line of vortices occur at ±45°. If this is not 
the case, and if pinning by twin boundaries is significant, 
then we should expect that the vortex lattice will align 
with the ionic lattice far from twin boundaries but will be 
deformed in the vicinity of a twin boundary in an effort to 
align itself with the twin boundary. On the other hand. 
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for e = 0, the vortex lattice would remain aligned with 
the twin boundaries everywhere except within vortex lat- 
tice domain boundaries which necessarily exist roughly 
midway between the twin boundaries. 

Neutron scattering experiments on YBCO [|| sug- 
gest that the vortex lattice is well-aligned with the twin 
boundaries and is close to being centered rectangular (the 
ratio of lattice constants is about 1.04) with /? « 73", 
with weak dependence on H . This corresponds to a ro- 
tation away from alignment with the ionic lattice by 9°. 
Four different orientational domains, related by reflec- 
tion in the (1, 1,0) axis and 90° rotation were reported. 
STM imaging of the YBCO vortex lattice also suggests 
that the (highly disordered) lattice has approximately 
centered rectangular symmetry with (3 ~ 77°. However, 
no evidence for the 9° tilt into alignment with the twin 
boundaries was reported. 

These neutron scattering results can be rather well fit- 
ted by the basic London model (e = 0) with mass 
anisotropy. There is a unique stretched hexagonal lattice 
which is aligned with either the (1, 1, 0) or (1, 1, 0) twin 
boundaries. For \x/\ = 1.5, a value roughly consistent 
with infrared and microwave experiments, this lattice has 
about the right shape. Taking into account the two crys- 
tallographic domains (related by interchanging Ax and 
Aj,) there are all together four vortex lattice domains, 
as seen experimentally. The experimental fact that the 
vortex lattice appears to be well aligned with the twin 
boundaries suggests that the tendency to align with the 
ionic lattice is small. No evidence for a bending of the 
vortex lattice (by 9°) into alignment with the ionic lattice 
far from the twin boundaries has so far been found. 

Low field Bitter decoration data on YBCO |^l[] show 
vortex lattice geometry with a very small distortion from 
hexagonal, consistent with a much smaller anisotropy 
Ax/Ay = 1.11 — 1.15. One may be tempted to attribute 
this apparent field dependence of /3 to the effects dis- 
cussed above in connection with Fig. |^(b). An alterna- 
tive, and perhaps more likely explanation, is a poor qual- 
ity of samples used in the Bitter decoration experiments 
that may have resulted in partial washing out of the a-b 
anisotropy otherwise present in clean crystals |^ . 

^SR experiments measure the field distribution 
P{B) = {l/VL) J d[B - B{f)]dxdy. This is shown in Fig. 
|l|(b). For (i ^ 60°, B{r) has two inequivalent saddle 
points leading to two peaks in P{B). P{B) is unaf- 
fected by effective mass anisotropy, as can be shown by 
the rescaling transformation, mentioned above. Existing 
liSR experiments show only a single peak p3| . 

The weak field dependence of (3, the alignment with 
twin boundaries and the single peak in P{B) all suggest 
that e is small in YBCO and that the normal London 
model, together with twin boundary pinning, provides a 
good fit to the data. Bitter decoration data do not fit well 
into this overall picture and further experimental work, 
preferably on untwinned YBCO or other tetragonal su- 
perconductors, will probably be necessary to clarify the 
importance of square lattice anisotropy in high- Tc super- 



conductors. 

The general approach to vortex lattices introduced 
here may be extended to low temperatures, but then the 
free energy takes a quite different form which is non- 
analytic in B at T ^ [|l|,0- 
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